compute_log_ccps <- function(xj, x0, hid, year, f_x, save_W) {
  require(Matrix)
  require(grf)
  newx <- rbind(xj, x0)

  colnames(newx) <- state.vars

  ND <- predict(rrf, newdata = newx) #
  ND[ND < 0.0001 | is.na(ND)] <- 0.0001
  #
  log_ccps <- log(ND)
  log_ccps <- cbind(log_ccps, log_ccps)

  return( as.matrix(log_ccps) )
}

compute_log_ccps_sim <- function(xj, x0) {
  require(Matrix)
  newx <- cbind(rbind(xj, x0))
  colnames(newx) <- state.vars
  newx <- data.frame(newx)
  newx[,"SN_ROR"] <- "917"
  newx[,"year"] <- "2009"

  newx[,"rent"] <- newx[,"rent0"]

  ccps <- predict.glm(ccp_model, newdata=newx, type="response")
  ccps[ccps < .0001] <- .0001 # 0.0001
  log_ccps <- log(ccps)
  return(as.matrix(log_ccps))
}


pre_compute_v_s <- function(
  x, # data for one household, all periods (excluding the last period, where the choice is not yet observed)
  state.vars, # names of state variables
  beta = .95, # discount factor
  nvars_u_c = 5, # number of common parameters in the flow utility function
  nvars_u_s = 4, # number of type-specific parameters in the flow utility function
  nvars_v = 8, # number of parameters in the terminal utility function
  sizes_sqm = c(39, 60, 76, 102), # unit size per # of rooms
  compute_psi = T,
  S, # number of unobserved states
  E = 50, # number of error pairs to be drawn
  save_W = T,
  base_path = "model_output",
  version=0,
  seed = 381271,
  g = 0
) {
  set.seed(seed)
  # Stored data transitions
  inc_transition_coefs <- as.numeric(readRDS(sprintf("%s/income_transition_coefs.RDS",base_path)))
  inc_transition_errors <- as.numeric(readRDS(sprintf("%s/income_transition_resid.RDS",base_path)))

  sav_transition_coefs <- as.numeric(readRDS(sprintf("%s/savings_transition_coefs.RDS",base_path)))
  sav_transition_errors <- as.numeric(readRDS(sprintf("%s/savings_transition_resid.RDS",base_path)))

  adu_transition_coefs <- as.matrix(readRDS(sprintf("%s/couple_transition_coefs.RDS",base_path)))

  kid_transition_coefs <- as.matrix(readRDS(sprintf("%s/children_transition_coefs.RDS",base_path)))


  start.time = Sys.time()
  require(tictoc)
  J <- 43 #  max(x$choice) # number of alternatives - 1
  x <- x[! is.na(x$choice),]

  hids <- unique(x$hid)
  all_t <- nrow(x)
  npar <- nvars_u_c + nvars_u_s + nvars_v * 3 + ifelse(S == 1,0,3)
  v <- matrix(data=NA, nrow=all_t * J * S, ncol=npar+6)
  t_index = 0


  rent_hs_default <- matrix(1, 1, 1)

  t_index <- 0
  tic()
  for(hid in hids) {

    data_n <- x[x$hid == hid,]
    years <- sort(data_n$year)
    max.t <- nrow(data_n)
    for(t in 1:max.t) {

      state <- data_n[data_n$year == years[t],state.vars]
      choice <- data_n[data_n$year == years[t],"choice"]

      # error draw
      ind.err <- sample(
        1:length(sav_transition_errors), E, replace=T
      )

      r <- mapped_v(
        as.numeric(state),
        beta,
        J,
        J-2, # terminal choice for period t+1
        nvars_u_c,
        nvars_u_s,
        nvars_v,
        compute_psi,
        inc_transition_coefs, # vector
        inc_transition_errors[ind.err], # vector
        sav_transition_coefs, # vector
        sav_transition_errors[ind.err], # vector
        adu_transition_coefs, # matrix
        kid_transition_coefs, # matrix
        hid,
        years[t],
        save_W,
        F, # simulation,
        1,
        S,
        sizes_sqm,
        matrix(1, 1, 1), # htype_hs default value
        matrix(0,4,10), # htype_hs_new default value
        g
      )
      v[1:(J*S) + t_index*(J*S),5:ncol(v)] <- r
      v[1:(J*S) + t_index*(J*S),1] <- hid
      v[1:(J*S) + t_index*(J*S),2] <- years[t]
      v[1:(J*S) + t_index*(J*S),3] <- choice
      v[1:(J*S) + t_index*(J*S),4] <- rep(1:J, S)

      t_index = t_index + 1
      if(t_index %% 10 == 0) cat(".")
      if(t_index %% 100 == 0) cat(" ")
      if(t_index %% 500 == 0) {
        pct_completed = 100*t_index/all_t
        time_elapsed <- as.numeric(as.difftime(Sys.time() - start.time), units="mins")
        time_remaining <- (time_elapsed / pct_completed) * (100 - pct_completed)
        cat(sprintf(" %.1f%%, ca. %.1f m remaining.\n", pct_completed, time_remaining))
      }
    }
  }
  cat(sprintf(" %d/%d observations\n",t_index,all_t))
  toc()
  cat("Done.")
  return(v)
}


#
update_q <- function(
  theta, # to calculate the log likelihood
  v, # hid, q and s are columns of vx
  full_data,
  state.vars,
  pi,
  q,
  base_path = base_path
) {

  J <- max(v[,"depvar"])
  hids <- unique(full_data$hid)
  q_new <- matrix(data=NA, nrow=length(hids), ncol=length(pi)+1)
  colnames(q_new) <- c("hid", sprintf("state%02d",1:length(pi)))
  q_new[,"hid"] <- hids

  expV <- exp(v[,names(theta)] %*% theta + v[,"psi"])
  expV[expV < 1e-16] <- 1e-16
  #i <- 0
  for(i in 1:length(hids)) {
    hid <- q_new[i,"hid"]
    choices_n <- full_data[full_data$hid == hid,"depvar"]
    choices_n <- as.numeric(choices_n[-1])

    exp_v <- expV[v[,"hid"] == hid]
    qn_new <- update_qn(
      choices_n, # sequence of choices for HH n, t = 1,..,T
      exp_v, # exp(v_jt)'s
      pi, # distribution of unobserved types in the data
      J # number of choices (ex. baseline choice, j=0)
    )
    if(any(is.nan(qn_new))) qn_new <- rep(1/8,8)
    q_new[i,sprintf("state%02d",1:length(pi))] <- qn_new
  }
  q_out <- q[,c("hid","year")]
  q_out <- merge(q_out, q_new, by="hid")
  return(q_out)
}


update_pi <- function(
  q
) {
  q.hh <- aggregate(
    q[,sprintf("state%02d", 1:(ncol(q)-2))],
    by = list(hid = q[,"hid"]),
    FUN=mean
  )
  q.hh <- q[,c("hid", sprintf("state%02d", 1:(ncol(q)-2)))]
  hids <- unique(q[,"hid"])
  pi_out <- apply(q.hh[2:ncol(q.hh)], 2, mean)
  return(pi_out)
}


get_LL <- function(
  v,
  theta,
  q0,
  start_index,
  J
) {
  numtheta <- as.numeric(theta)
  LL = LogLik_contribution(
    v,
    numtheta,
    as.matrix(q0[,2:ncol(q0)]),
    start_index,
    J
  )
  return(LL);
}


get_gradient <- function(
  v,
  theta,
  q0,
  start_index,
  J
) {
  numtheta <- as.numeric(theta)
  grad = gradient_contribution(
    v,
    numtheta,
    as.matrix(q0[,2:ncol(q0)]),
    start_index,
    J
  )
  return(grad);
}

update_theta <- function(
  v,
  data,
  theta,
  q0,
  J = 43,
  start_index,
  fixed
) {
  if(missing(fixed)) {
    loose <- theta
    fixed <- NULL
  } else {
    loose <- theta[!names(theta) %in% names(fixed)]
  }

  objective <- function(par) {
    this.theta <- c(par, fixed)
    this.theta <- this.theta[names(theta)]
    get_LL(v, this.theta, q0, start_index, J)
  }
  gradient <- function(par) {
    this.theta <- c(par, fixed)
    this.theta <- this.theta[names(theta)]
    grad <- get_gradient(v, this.theta, q0, start_index, J)
    grad[!names(theta) %in% names(fixed)]
  }
  fit <- optim(
    par = loose, fn = objective, gr = gradient, method = "BFGS",
    control = list(fnscale = -1, trace = 1, REPORT = 20, maxit=1000)
  )
  theta.fit <- fit$par
  names(theta.fit) <- names(loose)
  theta.fit <- c(theta.fit, fixed)
  theta.fit <- theta.fit[names(theta)]

  return(
    list(
      theta=theta.fit,
      LogLik = fit$value,
      convergence = fit$convergence,
      message = fit$message,
      counts = fit$counts,
      fixed = fixed
    )
  )
}



coeftest_boot <- function(
  result, # results object w/ data, v, theta
  seed=1029183, # seed
  B = 100, # replications
  E = 50,
  ntrees = 250,
  b.start = 1,
  load = F
) {
  base_path <- result$base_path
  theta <- result$theta
  state.vars <- result$state.vars

  nvars_u_c <- result$nvars_u_c
  nvars_u_s <- result$nvars_u_s
  nvars_v <- result$nvars_v

  v <- result$v
  v <- v[order(v[,"hid"], v[,"year"], v[,"s"], v[,"depvar"]),]


  if("pi" %in% names(result)) {
    S <- length(result$pi)
    q <- result$q
    q0 <- q[! duplicated(q[,1]),c(1,3:ncol(q))]
    q0 <- q0[order(q0[,"hid"]),]
    tmpfile <- "_het"
  } else  {
    S <- 1
    tmpfile <- ""
  }

  sizes_sqm = result$sizes_sqm
  if("g" %in% names(result)) g = result$g else g = 0

  J <- max(result$v[,"depvar"])
  if(load) { # load stored estimates only
    save_theta <- matrix(data=NA, nrow=B, ncol=length(theta))
    for(b in 1:B) {
      if(! file.exists(sprintf("%s/VCOV/bootstrap_tmp_b%03d.RDS", base_path, b))) {
        break;
        cat(sprintf("Bootstrap run %d missing. Please restart with b.start = %d.\n", b, b))
      } else {
        save_theta[b,]<- as.numeric(readRDS(sprintf("%s/VCOV/bootstrap_tmp_b%03d.RDS", base_path, b)))
      }
    }
    save_theta <- rbind(save_theta, theta)
  } else { # run bootstrap
    require(tictoc)
    require(grf)

    full_data <- result$data
    data <- full_data[ ! is.na(full_data$depvar_adj),]
    data$choice <- data$depvar_adj
    data <- data[order(data$hid, data$year),]

    beta <- result$beta
    sizes_sqm <- result$sizes_sqm


    hids <- unique(data$hid)
    N <- length(hids)

    save_theta <- matrix(data=NA, nrow=B, ncol=length(theta))

    for(b in b.start:B) {
      set.seed(seed + 2182 * b)
      cat("Repetition ", b, " - ", sep="")
      tic()

      draw <- sort(sample(1:N, N, replace=T))

      b.hids <- hids[draw]
      expand_draw <- unlist(lapply(
        b.hids, function(x) which(data$hid == x)
      ))

      draw <- data.frame(table(draw))
      names(draw)[2] <- "w"
      draw$hid <- hids[draw$draw]

      data_b <- data[data$hid %in% draw$hid,]


      #
      # # Subset of v
      vx_b <- v[v[,"hid"] %in% draw$hid,]

      if(S == 1) {
        q0_b <- draw[,c("hid", "w")] # use q0 as weights in update_theta()
        # # train random forest for choice J+1
        data_expand <- data[expand_draw,]
        choice <- data_expand$depvar_adj

        X <- data_expand[,state.vars]
        Y = as.integer(choice ==  (J - 2))
        rrf <<- regression_forest(
          X=X,
          Y=Y,
          num.trees = ntrees,
          tune.parameters="all"
        )
        print(rrf)
        #
        # # Compute v & psi
        vx_b <- pre_compute_v_s(
          x = data_b, # data for one household, all periods (excluding the last period, where the choice is not yet observed)
          state.vars = state.vars, # names of state variables
          beta = beta, # discount factor
          nvars_u_c = nvars_u_c, # number of constant variables in the flow utility function
          nvars_u_s = nvars_u_s, # number of type-specific variables in the flow utility function
          nvars_v = nvars_v, # number of variables in the terminal utility functions
          sizes_sqm = sizes_sqm,
          compute_psi = T, # set to F if just changing the utility function (skips psi function evaluation)
          S = S, # number of unobserved states
          E = E,
          save_W = F, # deprecated
          base_path = base_path,
          version = version,
          seed = seed,
          g = g
        )
        colnames(vx_b) <- colnames(v)
      } else {

        q0_b <- q0[q0[,"hid"] %in% draw$hid,]
        draw <- draw[order(draw[,"hid"]),]

        for(s in 1:S) q0_b[,sprintf("state%02d",s)] <- q0_b[,sprintf("state%02d",s)] * draw[,"w"]
        q0_b <- q0_b[,c("hid", sprintf(sprintf("state%02d",1:S)))]

      }
      # Get theta
      start_index <- which( ! duplicated(vx_b[,"hid"])) - 1
      start_index <- c(start_index, nrow(vx_b))
      fit <- update_theta(
        v = vx_b,
        data = data_b,
        theta = theta,
        q0 = q0_b,
        J = J,
        start_index = start_index,
        fixed = NULL
      )

      # Store estimates
      saveRDS(as.numeric(fit$theta), sprintf("%s/VCOV/bootstrap_tmp_b%03d.RDS", base_path, b))
      cat("\n")
      toc()
      cat("\n\n")
    }
    for(b in 1:B) {
      if(! file.exists(sprintf("%s/VCOV/bootstrap_tmp_b%03d.RDS", base_path, b))) {
        break;
        cat(sprintf("Bootstrap run %d missing. Please restart with b.start = %d.\n", b, b))
      } else {
        save_theta[b,]<- as.numeric(readRDS(sprintf("%s/VCOV/bootstrap_tmp_b%03d.RDS", base_path, b)))
      }
    }
    save_theta <- rbind(save_theta, theta)
  }

  df <- data.frame(
    Coef = theta,
    SE = apply(save_theta, 2, sd)
  )
  df$zval = df$Coef / df$SE
  df$pval = pnorm( -abs(df$zval) ) * 2
  return(
    list(
      CoefTable = df,
      theta_b = save_theta,
      seed=seed
    )
  )
}

# Simulation ####

update_rent0 <- function(htype, size, rent_htype_s, sizes_sqm, stay_length=0, growth_rate=0) {
  rent_htype_s[cbind(size, htype)] * sizes_sqm[size] * (1+growth_rate)^(-stay_length)
}



get_ccps <- function(
  new_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  rent_htype_s_new,
  compute_psi = T,
  i = NULL,
  supply_shock = NULL,
  pi_s = NULL,
  g = 0
) {


  npar <- nvars_u_c + nvars_u_s + 3*nvars_v + ifelse(is.null(pi_s),0,3)

  if(! is.null(pi_s)) {
    pi_s <- as.numeric(new_data[,sprintf("state%02d", 1:length(pi_s))])
  }


  r <- mapped_v(
    as.numeric(new_data[,state.vars]),
    beta,
    J,
    J-2, # terminal choice for period t+1
    nvars_u_c,
    nvars_u_s,
    nvars_v,
    compute_psi, # compute psi
    inc_transition_coefs, # vector
    inc_transition_errors, # vector
    sav_transition_coefs, # vector
    sav_transition_errors, # vector
    adu_transition_coefs, # matrix
    kid_transition_coefs, # matrix
    new_data[,"hid"],
    new_data[,"year"],
    F, # save_W
    T, # simulation,
    1,
    ifelse(is.null(pi_s),1,length(pi_s)),
    sizes_sqm,
    rent_htype_s,
    rent_htype_s_new,
    g
  )

  if(is.null(pi_s)) {
    v <- r[,1:npar] %*% as.numeric(theta) + r[,npar+1]
    if(! is.null(supply_shock)) {
      v[42] <- v[42] - supply_shock
    }
    v <- c(1, exp(v))


    ccp <- v/sum(v)
  } else {
    ccp <- 0
    for(s in 1:length(pi_s)) {

      v <- r[1:J + J*(s-1),1:npar] %*% theta + r[1:J,npar+1] #
      if(! is.null(supply_shock)) {
        v[42] <- v[42] - supply_shock
      }
      v <- c(1, exp(v))

      ccp <- pi_s[s] * v/sum(v) + ccp
    }
  }
  return(ccp)
}




get_d_ccps <- function(
  new_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors_sel, # vector
  sav_transition_coefs, # vector
  sav_transition_errors_sel, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  p42_coefs, # vector
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  rent_htype_s_new,
  compute_psi = T,
  i = NULL,
  supply_shock = NULL,
  pn, # ccps
  g = 0
) {
  npar <- nvars_u_c + nvars_u_s + 3*nvars_v
  d_pn <- dv_ccp(
    as.numeric(new_data[,state.vars]),
    beta,
    J,
    J-2, # terminal choice for period t+1
    nvars_u_c,
    nvars_u_s,
    nvars_v,
    compute_psi, # compute psi
    inc_transition_coefs, # vector
    inc_transition_errors_sel, # vector
    sav_transition_coefs, # vector
    sav_transition_errors_sel, # vector
    adu_transition_coefs, # matrix
    kid_transition_coefs, # matrix
    p42_coefs, # vector
    new_data[,"hid"],
    new_data[,"year"],
    F, # save_W
    T, # simulation,
    1,
    1,
    sizes_sqm,
    rent_htype_s,
    rent_htype_s_new,
    as.numeric(theta),
    pn,
    g
  )
  return(d_pn)
}


equil_objective <- function(
  base_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  rent_htype_s_new,
  compute_psi=F,
  ccps_seed,
  weights,
  int_demand_adj_fac,
  w_long_move0,
  buy_new_supply,
  new_supply_units,
  renov_prob_mat,
  rent_v,
  unit_sizes,
  supply_shock,
  supply_existing_baseline,
  supply_new_baseline,
  pi_s = NULL,
  target = 1e-08,
  g = 0
) {
  ccps <- matrix(NA, nrow=nrow(base_data), ncol=J+1)
  set.seed(ccps_seed)
  for(i in 1:nrow(base_data)) {
    ind.err <- sample(
      1:length(sav_transition_errors), E, replace=T
    )

    ccps[i,] <- get_ccps(
      new_data = base_data[i,],
      theta = theta,
      nvars_u_c = nvars_u_c,
      nvars_u_s = nvars_u_s,
      nvars_v = nvars_v,
      base_path = base_path,
      J = J,
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      state.vars = state.vars,
      sizes_sqm = sizes_sqm,
      beta = beta,
      rent_htype_s = rent_htype_s,
      rent_htype_s_new = rent_htype_s_new,
      compute_psi = compute_psi, #
      i = i,
      supply_shock = supply_shock,
      pi_s = pi_s,
      g = g
    )
  }
  ccps_save <- ccps

  #
  #
  for(j in 1:ncol(ccps)) ccps[,j] <- ccps_save[,j] * weights
  #

  # Fix supply of existing units to 'supply_existing' and re-scale other choices
  if(! is.null(supply_existing_baseline)) {
    supply_existing <- sum(ccps[,J-1])
    adj_sup <- (supply_existing_baseline/supply_existing) # scale proportionally
    ccps[,J-1] <- ccps[,J-1] * adj_sup

    tot_oth <- apply(ccps[,c(1:(J-2),J,J+1)], 1, sum) # total other than scaled
    adj_ccp <- (weights - ccps[,J-1]) / tot_oth
    for(j in c(1:(J-2),J,J+1)) ccps[,j] <- ccps[,j] * adj_ccp
    #
    # Fix supply of existing units to 'supply_existing' and re-scale other choices
    if(! is.null(supply_new_baseline)) {
      supply_new <- sum(ccps[,J])
      adj_sup <- (supply_new_baseline/supply_new) # scale proportionally
      ccps[,J] <- ccps[,J] * adj_sup

      tot_oth <- apply(ccps[,c(1:(J-2),J+1)], 1, sum) # total other than scaled and existing
      adj_ccp <- (weights - ccps[,J] - ccps[,J-1]) / tot_oth
      for(j in c(1:(J-2),J+1)) ccps[,j] <- ccps[,j] * adj_ccp
      #
    }
  }

  # Supply constraint: cannot accommodate more than buy_new_supply
  shock_size <- sum(ccps[,J]) / buy_new_supply - 1
  cat("Shock size =", round(shock_size, 6), "\n")
  #
  int_demand <- apply(ccps[,2:(J-2)], 2, sum)
  #
  # Inflow
  # w_long_move <- w_long_move0 # z_n
  ext_ccps <- ccps_save[,2:(J-2)]
  #
  #
  ext_demand <- t(ext_ccps) %*% w_long_move0

  # Total demand
  D <- int_demand_adj_fac * int_demand + ext_demand

  # Supply side: How many of the vacant units are upgraded to higher quality?
  # Vacancies
  vac_prob <- apply(ccps[,2:(J+1)],1,sum)

  vacant <- t(renov_prob_mat) %*% vac_prob
  #
  # Total supply
  S <- vacant + new_supply_units

  # Equilibrium
  diff <- D - S
  #
  Q <- sum(diff^2)
  print(Q)
  if(Q < target) {
    saveRDS(
      list(rent_fit = rent_htype_s, shock_size = shock_size),
      file=sprintf("%s/tmp_rent_fit.RDS", base_path)
    )
    Q <- 0
  }
  return(Q)
}





equil_grad <- function(
  base_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  p42_coefs, # vector
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  rent_htype_s_new,
  compute_psi=F, #
  ccps_seed,
  weights,
  int_demand_adj_fac,
  w_long_move0,
  buy_new_supply,
  new_supply_units,
  renov_prob_mat,
  rent_v,
  unit_sizes,
  supply_shock,
  supply_existing_baseline,
  supply_new_baseline,
  pi_s = NULL,
  g = 0
) {
  ccps <- matrix(NA, nrow=nrow(base_data), ncol=J+1)
  d_ccps <- array(NA, dim=c(nrow(base_data), J+1, 40)) #
  #
  set.seed(ccps_seed)
  for(i in 1:nrow(base_data)) {
    ind.err <- sample(
      1:length(sav_transition_errors), E, replace=T
    )
    ccps[i,] <- get_ccps(
      new_data = base_data[i,],
      theta = theta,
      nvars_u_c = nvars_u_c,
      nvars_u_s = nvars_u_s,
      nvars_v = nvars_v,
      base_path = base_path,
      J = J,
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      state.vars = state.vars,
      sizes_sqm = sizes_sqm,
      beta = beta,
      rent_htype_s = rent_htype_s,
      rent_htype_s_new = rent_htype_s_new,
      compute_psi = compute_psi, # update psi
      i = i,
      supply_shock = supply_shock,
      pi_s = pi_s,
      g = g
    )
    d_ccps[i,,] <- get_d_ccps(  #
      new_data = base_data[i,],
      theta = theta,
      nvars_u_c = nvars_u_c,
      nvars_u_s = nvars_u_s,
      nvars_v = nvars_v,
      base_path = base_path,
      J = J,
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      p42_coefs, # vector
      state.vars = state.vars,
      sizes_sqm = sizes_sqm,
      beta = beta,
      rent_htype_s = rent_htype_s,
      rent_htype_s_new = rent_htype_s_new,
      compute_psi = compute_psi, # update psi
      i = i,
      supply_shock = supply_shock,
      pn = ccps[i,],
      g = g
    )
  }
  ccps_save <- ccps
  for(j in 1:ncol(ccps)) ccps[,j] <- ccps_save[,j] * weights
  #
  # Fix supply of existing units to 'supply_existing' and re-scale other choices
  if(! is.null(supply_existing_baseline)) {
    other_js <- c(1:(J-2),J,J+1)
    #

    supply_existing <- sum(ccps[,J-1])
    d_supply_existing <- apply(d_ccps[,J-1,], 2, sum)

    adj_sup <- (supply_existing_baseline/supply_existing) # scale proportionally
    d_adj_sup <- -supply_existing_baseline * d_supply_existing / (supply_existing^2)

    ccps[,J-1] <- ccps[,J-1] * adj_sup

    tot_oth <- apply(ccps[,other_js], 1, sum) # total other than scaled
    adj_ccp <- (weights - ccps[,J-1]) / tot_oth
    for(j in other_js) ccps[,j] <- ccps[,j] * adj_ccp
    #
    d_adj_ccp <- d_tot_oth <- matrix(0, nrow=nrow(ccps), ncol=dim(d_ccps)[3])
    for(j in other_js) {
      d_tot_oth <- d_tot_oth + d_ccps[,j,]
    }
    denominator <- tot_oth^2
    for(k in 1:dim(d_ccps)[3]) {
      numerator1 <- (d_adj_sup[k] * ccps[,J-1] + adj_sup * d_ccps[,J-1,k]) * tot_oth
      numerator2 <- (weights - ccps[,J-1]) * d_tot_oth[,k]
      d_adj_ccp[,k] <- -(numerator1 + numerator2)/denominator
    }

    for(j in other_js) for(k in 1:dim(d_ccps)[3]) {
      d_ccps[,j,k] <- adj_ccp * d_ccps[,j,k] + d_adj_ccp[,k] * ccps[,j]
    }
    for(k in 1:dim(d_ccps)[3]) {
      d_ccps[,J-1,k] <- d_adj_sup[k] * ccps[,J-1] + adj_sup * d_ccps[,J-1,k]
    }


    if(! is.null(supply_new_baseline)) {
      other_js <- c(1:(J-2),J+1)
      #

      supply_new <- sum(ccps[,J])
      d_supply_new <- apply(d_ccps[,J,], 2, sum)

      adj_sup <- (supply_new_baseline/supply_new) # scale proportionally
      d_adj_sup <- -supply_new_baseline * d_supply_new / (supply_new^2)

      ccps[,J] <- ccps[,J] * adj_sup

      tot_oth <- apply(ccps[,other_js], 1, sum) # total other than scaled
      adj_ccp <- (weights - ccps[,J] - ccps[,J-1]) / tot_oth
      for(j in other_js) ccps[,j] <- ccps[,j] * adj_ccp
      #
      d_adj_ccp <- d_tot_oth <- matrix(0, nrow=nrow(ccps), ncol=dim(d_ccps)[3])
      for(j in other_js) {
        d_tot_oth <- d_tot_oth + d_ccps[,j,]
      }
      denominator <- tot_oth^2
      for(k in 1:dim(d_ccps)[3]) {
        numerator1 <- (d_adj_sup[k] * ccps[,J] + adj_sup * d_ccps[,J,k]) * tot_oth
        numerator2 <- (weights - ccps[,J] - ccps[,J-1]) * d_tot_oth[,k]
        d_adj_ccp[,k] <- -(numerator1 + numerator2)/denominator
      }

      for(j in other_js) for(k in 1:dim(d_ccps)[3]) {
        d_ccps[,j,k] <- adj_ccp * d_ccps[,j,k] + d_adj_ccp[,k] * ccps[,j]
      }
      for(k in 1:dim(d_ccps)[3]) {
        d_ccps[,J,k] <- d_adj_sup[k] * ccps[,J] + adj_sup * d_ccps[,J,k]
      }
    }
  }
  #
  # 1. S(r) - D(r)
  #
  # 1.1 S(r)
  # Vacancies
  vac_prob <- apply(ccps[,2:(J+1)],1,sum)
  vacant <- t(renov_prob_mat) %*% vac_prob
  #
  # Total supply
  S <- vacant + new_supply_units

  # 1.2 D(r)
  # Internal
  int_demand <- apply(ccps[,2:(J-2)], 2, sum)
  # Inflow
  ext_ccps <- ccps_save[,2:(J-2)]
  ext_demand <- t(ext_ccps) %*% w_long_move0

  # Total demand
  D <- int_demand_adj_fac * int_demand + ext_demand

  # 2. dS(r) / dr - dD(r) / dr
  dS_dr <- dD_dr <- matrix(NA, length(S), length(S))

  # Jacobiab:  df_j / dx_k = J(j,k)

  # dS_j(r) / dr_k = dS_sec_j(r) / dr_k
  for(j in 1:length(S)) for(k in 1:length(S)) {
      component = weights * renov_prob_mat[,j] * d_ccps[,1,k] # check dims / assignments!
      dS_dr[j,k] = -sum(component) # check dims / assignments!
  }
  # 2.2 dD_j(r) / dr_k
  for(j in 1:length(S)) for(k in 1:length(S)) {
    component = (int_demand_adj_fac[j] * weights + w_long_move0)  * d_ccps[,j+1,k] # check dims / assignments!
    dD_dr[j,k] = sum(component) # check dims / assignments!
  }

  # 3. dQ / dr
  diff <- S - D
  d_diff <- dS_dr - dD_dr
  Q <- sum(diff^2)
  dQ <- 2 * t(d_diff) %*% diff

  return(dQ)
}


equil_quant <- function(
  base_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  compute_psi=F,
  ccps_seed,
  weights,
  int_demand_adj_fac,
  w_long_move0,
  buy_new_supply,
  new_supply_units,
  renov_prob_mat,
  rent_v,
  unit_sizes,
  supply_shock,
  supply_existing_baseline,
  supply_new_baseline,
  pi_s = NULL,
  g = 0
) {
  rent_htype_s_new <- rent_htype_s
  ccps <- matrix(NA, nrow=nrow(base_data), ncol=J+1)
  set.seed(ccps_seed)
  for(i in 1:nrow(base_data)) {
    ind.err <- sample(
      1:length(sav_transition_errors), E, replace=T
    )

    ccps[i,] <- get_ccps(
      new_data = base_data[i,],
      theta = theta,
      nvars_u_c = nvars_u_c,
      nvars_u_s = nvars_u_s,
      nvars_v = nvars_v,
      base_path = base_path,
      J = J,
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      state.vars = state.vars,
      sizes_sqm = sizes_sqm,
      beta = beta,
      rent_htype_s = rent_htype_s,
      rent_htype_s_new = rent_htype_s_new,
      compute_psi = compute_psi, # update psi
      i = i,
      supply_shock = supply_shock,
      pi_s = pi_s,
      g = g
    )
  }
  ccps_save <- ccps
  #cat("\n")

  #
  #
  for(j in 1:ncol(ccps)) ccps[,j] <- ccps_save[,j] * weights
  #

  # Fix supply of existing units to 'supply_existing' and re-scale other choices
  if(! is.null(supply_existing_baseline)) {
    supply_existing <- sum(ccps[,J-1])
    adj_sup <- (supply_existing_baseline/supply_existing) # scale proportionally
    ccps[,J-1] <- ccps[,J-1] * adj_sup

    tot_oth <- apply(ccps[,c(1:(J-2),J,J+1)], 1, sum) # total other than scaled
    adj_ccp <- (weights - ccps[,J-1]) / tot_oth
    for(j in c(1:(J-2),J,J+1)) ccps[,j] <- ccps[,j] * adj_ccp
    #
    if(! is.null(supply_new_baseline)) {
      supply_new <- sum(ccps[,J])
      adj_sup <- (supply_new_baseline/supply_new) # scale proportionally
      ccps[,J] <- ccps[,J] * adj_sup

      tot_oth <- apply(ccps[,c(1:(J-2),J+1)], 1, sum) # total other than scaled and existing
      adj_ccp <- (weights - ccps[,J] - ccps[,J-1]) / tot_oth
      for(j in c(1:(J-2),J+1)) ccps[,j] <- ccps[,j] * adj_ccp
      #
    }
  }
  # Supply constraint: cannot accommodate more than buy_new_supply
  shock_size <- sum(ccps[,J]) / buy_new_supply - 1
  cat("Shock size =", round(shock_size, 6), "\n")
  #
  int_demand <- apply(ccps[,2:(J-2)], 2, sum)
  #
  # Inflow
  # w_long_move <- w_long_move0 # z_n
  ext_ccps <- ccps_save[,2:(J-2)]
  #
  #
  ext_demand <- t(ext_ccps) %*% w_long_move0

  # Total demand
  D <- int_demand_adj_fac * int_demand + ext_demand

  # Total supply
  vac_prob <- apply(ccps[,2:(J+1)],1,sum)

  vacant <- t(renov_prob_mat) %*% vac_prob
  #
  S <- vacant + new_supply_units

  return(list(D = D, S = S, ds = shock_size))
}



update_rent_psqm <- function(
  base_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  p42_coefs, # vector
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  ccps_seed,
  weights,
  int_demand_adj_fac,
  w_long_move0,
  buy_new_supply,
  new_supply_units,
  renov_prob_mat,
  unit_sizes,
  method = "BFGS",
  supply_shock = 1e-3,
  supply_existing_baseline,
  supply_new_baseline,
  pi_s = NULL,
  target=1e-08,
  max.iter = 100,
  g = 0
) {
  if(missing(supply_existing_baseline)) supply_existing_baseline <- NULL
  if(missing(supply_new_baseline)) supply_new_baseline <- NULL

  par_init <- as.vector(t(rent_htype_s))
  objective <- function(par) {
    rent_htype_s <- matrix(par, 4, 10, byrow = T)
    rent_htype_s_new <- rent_htype_s
    rent_v <- par * unit_sizes
    equil_objective(
      base_data=base_data,
      theta=theta,
      nvars_u_c=nvars_u_c,
      nvars_u_s=nvars_u_s,
      nvars_v=nvars_v,
      base_path=base_path,
      J=J,
      inc_transition_coefs=inc_transition_coefs, # vector
      inc_transition_errors=inc_transition_errors, # vector
      sav_transition_coefs=sav_transition_coefs, # vector
      sav_transition_errors=sav_transition_errors, # vector
      adu_transition_coefs=adu_transition_coefs, # matrix
      kid_transition_coefs=kid_transition_coefs, # matrix
      state.vars=state.vars,
      sizes_sqm=sizes_sqm,
      beta=beta,
      rent_htype_s=rent_htype_s,
      rent_htype_s_new=rent_htype_s_new,
      compute_psi=T, # update psi every 10 iterations only
      ccps_seed=ccps_seed,
      weights=weights,
      int_demand_adj_fac=int_demand_adj_fac,
      w_long_move0=w_long_move0,
      buy_new_supply=buy_new_supply,
      new_supply_units=new_supply_units,
      renov_prob_mat=renov_prob_mat,
      rent_v = rent_v,
      unit_sizes = unit_sizes,
      supply_existing_baseline = supply_existing_baseline,
      supply_new_baseline = supply_new_baseline,
      supply_shock = supply_shock,
      pi_s = pi_s,
      target = target,
      g = g
    ) * 100
  }
  gradient <- function(par) {
    rent_htype_s <- matrix(par, 4, 10, byrow = T)
    rent_htype_s_new <- rent_htype_s
    rent_v <- par * unit_sizes
    equil_grad(
      base_data=base_data,
      theta=theta,
      nvars_u_c=nvars_u_c,
      nvars_u_s=nvars_u_s,
      nvars_v=nvars_v,
      base_path=base_path,
      J=J,
      inc_transition_coefs=inc_transition_coefs, # vector
      inc_transition_errors=inc_transition_errors, # vector
      sav_transition_coefs=sav_transition_coefs, # vector
      sav_transition_errors=sav_transition_errors, # vector
      adu_transition_coefs=adu_transition_coefs, # matrix
      kid_transition_coefs=kid_transition_coefs, # matrix
      p42_coefs=p42_coefs, # vector
      state.vars=state.vars,
      sizes_sqm=sizes_sqm,
      beta=beta,
      rent_htype_s=rent_htype_s,
      rent_htype_s_new=rent_htype_s_new,
      compute_psi=T, # update psi every 10 iterations only
      ccps_seed=ccps_seed,
      weights=weights,
      int_demand_adj_fac=int_demand_adj_fac,
      w_long_move0=w_long_move0,
      buy_new_supply=buy_new_supply,
      new_supply_units=new_supply_units,
      renov_prob_mat=renov_prob_mat,
      rent_v = rent_v,
      unit_sizes = unit_sizes,
      supply_existing_baseline = supply_existing_baseline,
      supply_new_baseline = supply_new_baseline,
      supply_shock = supply_shock,
      pi_s = pi_s,
      g = g
    ) * 100
  }
  fit <- optim(
    par = par_init, fn = objective, gr = gradient, method = method,
    control = list(
      trace = 1, REPORT = 1, maxit=max.iter, abstol=1e-12
       #,parscale=rep(1, 40)
    )
  )
  rent_htype_s <- matrix(fit$par, 4, 10, byrow = T)

  return(
    list(
      rent_htype_s=rent_htype_s,
      LogLik = fit$value,
      convergence = fit$convergence,
      message = fit$message,
      counts = fit$counts
    )
  )
}



sim_welfare <- function(
  base_data,
  theta,
  nvars_u_c,
  nvars_u_s,
  nvars_v,
  base_path,
  J,
  inc_transition_coefs, # vector
  inc_transition_errors, # vector
  sav_transition_coefs, # vector
  sav_transition_errors, # vector
  adu_transition_coefs, # matrix
  kid_transition_coefs, # matrix
  state.vars,
  sizes_sqm,
  beta,
  rent_htype_s,
  rent_htype_s_post,
  compute_psi=T,
  ccps_seed,
  supply_shock,
  E = 50,
  n_draws = 50,
  n_draw_seed = 1319238,
  pi_s = NULL,
  g = 0
) {
  npar <- nvars_u_c + nvars_u_s + 3*nvars_v + ifelse(is.null(pi_s),0,3)
  rent_htype_s_new <- rent_htype_s
  rent_htype_s_new_post <- rent_htype_s_post


  if(is.null(pi_s)) {
    save_v <- save_v_post <- matrix(NA, nrow=nrow(base_data), ncol=J)
    save_v_hinc <- save_v_base <- matrix(NA, nrow=nrow(base_data), ncol=J+1)
  } else {
    save_v <- save_v_post <- array(NA, dim=c(nrow(base_data), length(pi_s), J))
    save_v_hinc <- save_v_base <- array(NA, dim=c(nrow(base_data), length(pi_s), J+1))
  }


  set.seed(ccps_seed)
  for(i in 1:nrow(base_data)) {
    ind.err <- sample(
      1:length(sav_transition_errors), E, replace=T
    )
    new_data <- base_data[i,state.vars]
    r <- mapped_v(
      as.numeric(new_data),
      beta,
      J,
      J-2, # terminal choice for period t+1
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      compute_psi, # compute psi
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      base_data[i,"hid"],
      base_data[i,"year"],
      F, # save_W
      T, # simulation,
      1,
      ifelse(is.null(pi_s),1,length(pi_s)),
      sizes_sqm,
      rent_htype_s,
      rent_htype_s_new,
      g
    )
    r_base_diff <- r
    r_base <- mapped_v_total( # v_j // mapped_v w/o subtracting v_0
      as.numeric(new_data),
      beta,
      J,
      J-2, # terminal choice for period t+1
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      compute_psi, # compute psi
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      base_data[i,"hid"],
      base_data[i,"year"],
      F, # save_W
      T, # simulation,
      1,
      ifelse(is.null(pi_s),1,length(pi_s)),
      sizes_sqm,
      rent_htype_s,
      rent_htype_s_new,
      g
    )

    r_post <- mapped_v(
      as.numeric(new_data),
      beta,
      J,
      J-2, # terminal choice for period t+1
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      compute_psi, # compute psi
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      base_data[i,"hid"],
      base_data[i,"year"],
      F, # save_W
      T, # simulation,
      1,
      ifelse(is.null(pi_s),1,length(pi_s)),
      sizes_sqm,
      rent_htype_s_post,
      rent_htype_s_new_post,
      g
    )
    new_data["income"] <- new_data["income"] + 1
    r_hinc_diff <- mapped_v( # diff: v_j - v_0
      as.numeric(new_data),
      beta,
      J,
      J-2, # terminal choice for period t+1
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      compute_psi, # compute psi
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      base_data[i,"hid"],
      base_data[i,"year"],
      F, # save_W
      T, # simulation,
      1,
      ifelse(is.null(pi_s),1,length(pi_s)),
      sizes_sqm,
      rent_htype_s,
      rent_htype_s_new,
      g
    )
    r_hinc <- mapped_v_total( # v_j // mapped_v w/o subtracting v_0
      as.numeric(new_data),
      beta,
      J,
      J-2, # terminal choice for period t+1
      nvars_u_c,
      nvars_u_s,
      nvars_v,
      compute_psi, # compute psi
      inc_transition_coefs, # vector
      inc_transition_errors[ind.err], # vector
      sav_transition_coefs, # vector
      sav_transition_errors[ind.err], # vector
      adu_transition_coefs, # matrix
      kid_transition_coefs, # matrix
      base_data[i,"hid"],
      base_data[i,"year"],
      F, # save_W
      T, # simulation,
      1,
      ifelse(is.null(pi_s),1,length(pi_s)),
      sizes_sqm,
      rent_htype_s,
      rent_htype_s_new,
      g
    )
    # (r_hinc - r_hinc_diff)*theta = v_j - (v_j - v_0) = v_0

    if(is.null(pi_s)) {
      save_v[i,] <- r[,1:npar] %*% as.numeric(theta) + r[,npar+1]
      save_v_post[i,] <- r_post[,1:npar] %*% as.numeric(theta) + r_post[,npar+1]

      save_v_hinc[i,2:(J+1)] <- r_hinc[,1:npar] %*% as.numeric(theta) + r_hinc[,npar+1]
      save_v_hinc[i,1] <- (r_hinc[1,1:npar] - r_hinc_diff[1,1:npar]) %*% as.numeric(theta) +
        (r_hinc[1,npar+1] - r_hinc_diff[1,npar+1])

      save_v_base[i,2:(J+1)] <- r_base[,1:npar] %*% as.numeric(theta) + r_base[,npar+1]
      save_v_base[i,1] <- (r_base[1,1:npar] - r_base_diff[1,1:npar]) %*% as.numeric(theta) +
        (r_base[1,npar+1] - r_base_diff[1,npar+1])

      if(! is.null(supply_shock)) {
        save_v_post[i,42] <- save_v_post[i,42] - supply_shock
      }
    } else {
      for(s in 1:length(pi_s)) {
        v <- r[1:J + J*(s-1),1:npar] %*% theta + r[1:J,npar+1] # (assuming psi does not depend on s, for simplicity)
        v_post <- r_post[1:J + J*(s-1),1:npar] %*% theta + r_post[1:J,npar+1]

        v_hinc <- r_hinc[1:J + J*(s-1),1:npar] %*% theta + r_hinc[1:J,npar+1]
        v_hinc <- c(
          (r_hinc[1 + J*(s-1),1:npar] - r_hinc_diff[1 + J*(s-1),1:npar]) %*% theta +
            (r_hinc[1,npar+1] - r_hinc_diff[1,npar+1]), # v_0
          v_hinc # v_j's, j = 1 ,..., J
        )

        v_base <- r_base[1:J + J*(s-1),1:npar] %*% theta + r_base[1:J,npar+1]
        v_base <- c(
          (r_base[1 + J*(s-1),1:npar] - r_base_diff[1 + J*(s-1),1:npar]) %*% theta +
            (r_base[1,npar+1] - r_base_diff[1,npar+1]),
          v_base
        )

        if(! is.null(supply_shock)) {
          v_post[42] <- v_post[42] - supply_shock
        }
        save_v[i,s,] <- v
        save_v_post[i,s,] <- v_post

        save_v_hinc[i,s,] <- v_hinc
        save_v_base[i,s,] <- v_base
      }
    }
  }

  require(evd)
  vdiff <- vdiff_neg <- matrix(NA, nrow=nrow(base_data), ncol=n_draws)
  vdiff_inc <- matrix(NA, nrow=nrow(base_data), ncol=n_draws)

  # Draw s types
  set.seed(n_draw_seed)
  if(length(pi_s)>1) {
    s_draw <- numeric(nrow(base_data))
    for(i in 1:nrow(base_data)) {
      this_prob <- as.numeric(base_data[i,sprintf("state%02d", 1:length(pi_s))])
      s_draw[i] <- sample(1:length(pi_s), 1, prob=this_prob)
    }
  }

  set.seed(n_draw_seed)
  for(ndraw in 1:n_draws) {
    error_draws <- matrix(
      data = rgumbel(nrow(base_data)*J),
      nrow=nrow(base_data),
      ncol=J
    )
    v <- v_post <- matrix(NA, nrow=nrow(base_data), ncol=J)
    v_hinc <- v_base <- matrix(NA, nrow=nrow(base_data), ncol=J+1)

    if(length(pi_s)>1) {
      for(i in 1:nrow(base_data)) {

        v[i,] <- save_v[i,s_draw[i],] + error_draws[i,]
        v_post[i,] <- save_v_post[i,s_draw[i],] + error_draws[i,]

        v_hinc[i,] <- save_v_hinc[i,s_draw[i],] + c(0, error_draws[i,])
        v_base[i,] <- save_v_base[i,s_draw[i],] + c(0, error_draws[i,])
      }
    } else {
      v <- save_v + error_draws
      v_post <- save_v_post + error_draws

      v_hinc <- save_v_hinc + c(0, error_draws)
      v_base <- save_v_base + c(0, error_draws)
    }

    max_v <- apply(cbind(v, 0), 1, max)

    max_v_post <- apply(cbind(v_post,0), 1, max)

    max_v_hinc <- apply(v_hinc, 1, max)
    max_v_base <- apply(v_base, 1, max)

    vdiff[,ndraw] <- max_v_post - max_v
    vdiff_neg[,ndraw] <- as.integer(max_v_post < max_v)
    vdiff_inc[,ndraw] <- max_v_hinc - max_v_base
  }
  return(
    data.frame(
      vdiff = apply(-vdiff, 1, mean), # expected change in valuation
      vdiff_pos = apply(vdiff_neg, 1, mean), # expectation: increase in valuation negative (from supply reduction)
      vdiff_inc = apply(vdiff_inc, 1, mean) # change in valuation from 1 EUR additional income
    )
  )
}

